function [modelSolution,maxRelChange] = solveModelPde(xGrid,sigmaI,xDrift,xVolatility,params,algParams)
%SOLVEMODELPDE Solves the model globally on the specified state grid 
%   Inputs:
%     - xGrid: grid for state variable x (in practice, x=\tilde{\sigma}^2,
%              but this is not assumed)
%     - sigmaI: value of \tilde{\sigma}(x) on the grid
%     - xDrift: (arithmetic) drift of state variable x on the grid
%     - xVolatility: (arithmetic) volatility of state variable x on the grid
%     - params: struct array containing model parameters
%     - algParams: struct array containing algorithm parameters (requires
%                  the parameters nSolutionSteps, dt, convCrit)
%   Outputs:
%     - modelSolution: struct array containing model solution functions
%     - maxRelChange: convergence metric

  % unpack params
  [sigma,chi,gamma,rho,phi,iota0,delta,a0,g0,adjBondGrowth0,alphaA,alphaBondGrowth] = ...
    aux.unpack(params,'sigma','chi','gamma','rho','phi','iota0','delta','a0','g0','adjBondGrowth0','alphaA','alphaBondGrowth');
  % unpack algParams
  [nSolutionSteps,dt,convCrit] = aux.unpack(algParams,'nSolutionSteps','dt','convCrit');
  
  % compute linear rules for policy and exogenous variables
  a = a0 - alphaA*(sigmaI-sigma);
  g = g0*ones(size(xGrid));
  adjBondGrowth = adjBondGrowth0 + alphaBondGrowth*(sigmaI-sigma);
  
  % prepare solution iterative loop
  differentiator = finiteDifferentiator(xGrid,3);

  theta = 0.5*ones(size(xGrid));
  logValue = zeros(size(xGrid)); %log(v)

  diffusion = 1/2*xVolatility.^2;
  diffusion([1,end]) = 0;
  advection = xDrift;
  zeroInhomogeneity = zeros(size(xGrid));

  % main iterative loop to solve PDEs
  for i=1:nSolutionSteps
    % compute volatilities of theta and value
    sigmaTheta = differentiator.d1(log(theta)).*xVolatility;
    sigmaValue = differentiator.d1(logValue).*xVolatility;
    % derived quantities
    consWealthRatio = rho*ones(size(xGrid));
    totalWealth = (1 + phi*(a-g-iota0))./(1-theta+phi*consWealthRatio);
    consCapRatio = consWealthRatio.*totalWealth;
    aggregateRiskPremium = sigmaValue.*sigmaTheta;
    idiosyncraticRiskPremium = gamma*chi^2*(1-theta).^2.*sigmaI.^2;
    muTheta = consWealthRatio + adjBondGrowth - aggregateRiskPremium - idiosyncraticRiskPremium;
    iota = iota0 + ((1-theta).*(a-g-iota0)-consWealthRatio)./(1-theta+phi*consWealthRatio);
    growthRate = investmentTechnology(iota,phi,iota0) - delta;
    % pde coefficients
    discountingTheta = -muTheta;
    discountingLogValue = -rho*ones(size(xGrid));
    inhomogeneityLogValue = -(gamma-1)*(rho*log(consCapRatio) + growthRate - idiosyncraticRiskPremium/2) + sigmaValue.^2/2;

    % update theta
    thetaLong = linearParabolicPDE_timeStep([xGrid(1)-1;xGrid;xGrid(end)+1],[0;theta;0],...
      dt,[NaN;diffusion;NaN],[NaN;advection;NaN],[NaN;discountingTheta;NaN],[NaN;zeroInhomogeneity;NaN],0,0);
    thetaNew = thetaLong(2:end-1);
    % update logValue
    logValueLong = linearParabolicPDE_timeStep([xGrid(1)-1;xGrid;xGrid(end)+1],[0;logValue;0],...
      dt,[NaN;diffusion;NaN],[NaN;advection;NaN],[NaN;discountingLogValue;NaN],[NaN;inhomogeneityLogValue;NaN],0,0);
    logValueNew = logValueLong(2:end-1);

    % compute convergence metric
    relChangeTheta = (thetaNew-theta)./(abs(thetaNew)+abs(theta))*2;
    relChangeLogValue = (logValueNew-logValue)./(abs(logValueNew)+abs(logValue))*2;
    maxRelChangeTheta = max(abs(relChangeTheta));
    maxRelChangeLogValue = max(abs(relChangeLogValue));
    maxRelChange = max(maxRelChangeTheta,maxRelChangeLogValue);

    % transfer variables
    theta = thetaNew;
    logValue = logValueNew;

    if maxRelChange < convCrit
      break;
    end
  end
  
  % compute other variables
  consWealthRatio = rho*ones(size(xGrid));
  totalWealth = (1 + phi*(a-g-iota0))./(1 - theta + phi*consWealthRatio);
  consCapRatio = rho*totalWealth;
  qK = (1-theta).*totalWealth;
  qB = theta.*totalWealth;
  iota = iota0 + ((1-theta).*(a-g-iota0)-rho)./(1-theta+phi*consWealthRatio);
  growthRate = investmentTechnology(iota,phi,iota0) - delta; %this is the growth rate of capital
  sigmaTotalWealth = differentiator.d1(log(totalWealth)).*xVolatility;
  sigmaTheta = differentiator.d1(log(theta)).*xVolatility;
  sigmaValue = differentiator.d1(logValue).*xVolatility;
  aggrPriceOfRiskTotWealthNumeraire = -sigmaValue;
  aggrPriceOfRisk = aggrPriceOfRiskTotWealthNumeraire + sigmaTotalWealth;
  equityPremium = aggrPriceOfRisk.*(-sigmaTheta./(1-theta));
  equityExcRetVol =  sigmaTheta./(1-theta);
  govSpending = g;
  
  % risk-free rate and similar objects
  dqRel = differentiator.d1(log(totalWealth));
  d2qRel = differentiator.d2(log(totalWealth)) + dqRel.^2;
  muTotalWealth = dqRel.*xDrift + d2qRel.*xVolatility.^2/2;
  growthRateCons = muTotalWealth + growthRate;
  riskFreeRate = consWealthRatio + growthRateCons - aggrPriceOfRisk.*sigmaTotalWealth - gamma*chi^2*(1-theta).^2.*sigmaI.^2;
  expectedReturnBonds = riskFreeRate + aggrPriceOfRisk.*(sigmaTotalWealth + sigmaTheta);
  
  % save outputs in a struct array
  modelSolution = aux.pack(sigmaI,a,govSpending,theta,logValue,totalWealth,qK,qB,iota,growthRate,consCapRatio,aggrPriceOfRisk,equityPremium,equityExcRetVol,adjBondGrowth,growthRateCons,riskFreeRate,expectedReturnBonds,sigmaTotalWealth,sigmaTheta,sigmaValue);
  
end

